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Abstract 

Quantum-mechanical analysis based on an exact sum rule is used to extract 
an semiclassical angle-dependent energy function for transition metal ions in 
biomolecules. The angular dependence is simple but different from exist- 
ing classical potentials. Comparison of predicted energies with a computer- 
generated database shows that the semiclassical energy function is remarkably 
accurate, and that its angular dependence is optimal. 
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Biomolecular modeling with classical potentials has become an increasingly important 
tool in problems such as the determination of protein structure and function, and the design 
of new molecules with desired properties. With the continuing availability of increasingly 
powerful computers, one can only expect this growth to continue. A major roadblock toward 
expanded use of modeling with classical potentials is the absence of sufficiently reliable force 
laws for transition metals in biomolecules. The ability to treat transition metals is impor- 
tant because the active sites of many proteins are defined by transition metals; also smaller 
biologically active molecules often have transition metals as crucial constituents. In fact, 
transition metal complexes have been proposed as crucial ingredients in the origin of life it- 
self |]J . Unlike the s-p constituents of biomolecules, which usually have unique, well-defined 
bonding configurations (such as sp 2 planar coordination), transition metals can adopt a 
broad range of asymmetric environments. This asymmetry is often important for the func- 
tioning of enzymes. Thus, for biomolecular modeling, one needs a "generic" potential which 
treats essentially all physically reasonable environments instead of perturbations relative to 
as single well-defined structure. A priori, one does not know the functional form of such a 
generic potential. The pair approximation, which ignores angular constraints, is applicable 
to simple metal ions, but not to transition metals. The transition-metal d-orbitals lead 
to complex angular forces which are manifested, for example, in the frequent occurrence 
of Cu 2+ and Ni 2+ ions in square-planar environments that are unexpected on the basis of 
pair interactions alone. In existing simulation codes based on classical potentials, the an- 
gular terms are usually either ignored @-f|], on the assumption that direct ligand-ligand 
interactions can take up most of the "slack" , or they are treated with simple assumed an- 
gular forms. The latter range from quadratic or higher order expansions about observed 
equilibrium bond angles [|5|-§] to more sophisticated expansions in trigonometric functions 
P^T2|. However, there has been no derivation of the angular form of classical potentials 
from quantum mechanics. 

In this Letter, I use quantum-mechanical analysis to derive an energy function for d- 
electrons based on the local environment in biomolecules. The energy function has a "semi- 
classical" form, in the sense that it is slightly more complex than a classical additive sum of 
ligand-ligand interactions, but is still straightforward to treat in molecular modeling codes. 
To test the energy function, I generate a large number of random transition metal environ- 
ments and evaluate their exact energies as a test set. The ci-electron energy is described 
with surprising precision. The accuracy is much better than that of commonly used func- 
tional forms, and significantly improves on that of additive energy functions. The angular 
dependence of an energy function obtained by fitting to the exact energies is very similar to 
that derived analytically. 

In biomolecules, transition metals are typically in a "coordination" bonding configura- 
tion. This differs from metallic bond in elemental transition metals in that the <i-states 
usually hybridize with ligand orbitals at lower energies, rather than other d-orbitals at the 
same energy. This leads to well-defined discrete charge states. The physics of coordination 
bonding is well described by the ligand field theory [T^] (LFT), which treats the d-shell in 
a transition-metal ligand complex by an effective d-d Hamiltonian: 

H d = h^ u \df,){d u \ . (1) 
Here \dp) and \d u ) are <i-basis orbitals on the transition- metal ion, and the contain the 
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effects of the ligands in a perturbative fashion: 



(2) 



where the \i) are orbitals on the ligands that hybridize with the ci-shell, and E d and Ei are the 
c?-shell and ligand-orbital energies, respectively. (The ligands are all taken to be equivalent 
for simplicity, but the more general cases are treated straightforwardly.) This approximate 
treatment describes the systematics of ci-shell splittings in transition metal complexes quite 
well, although the electronic transition energies are not obtained quantitatively. In the case 
where only cr-type interactions between the ligands and the <i-shell are present, the matrix 
elements of the effective Hamiltonian can be written [|16H as 
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where the Y u have the angular dependence of the <i-basis orbitals, and the radial function 
e(rj) includes the effects of the energy denominator as well as the matrix elements. 

The ci-electron energy associated with Hd is obtained by simply adding the eigenenergies 
of the occupied d-states. This approximation is justified when comparing structural energies 
within a single well-defined charge/spin state. We focus on the "ligand-field stabilization 
energy" -Elfse = J2 n e n — N d s. Here the first term denotes the eigenvalue sum, N d is 
the number of (^-electrons, and e is the average energy of the (i-complex (including both 
occupied and unoccupied states). As indicated in Fig. 1, splitting of the <i-complex by 
ligand-field interactions provides a negative (stabilizing) contribution to -Elfse if the <i-shell 
is partly filled. The stabilizing contribution is enhanced if there is a gap between the highest 
occupied and lowest unoccupied states, as occurs for Cu 2+ and Ni 2+ ions in the square-planar 
coordination. We define the half- width W of the <i-complex as the rms deviation of the energy 
eigenvalues from the (i-complex average energy e. In the first approximation, one expects 
that -Elfse should be proportional to W. 

The (i-electron energy function developed here gives -Elfse &s & 

simple function of the 

ligand positions. It is based on an exact sum rule that for W. Explicit calculation via 
Eqs. ([I]) and (|2|) shows that 



5W 2 = E(^n - ef = Tr(H d - el) 2 = £ U„ 



(4) 



h3 



where the ligand-ligand interaction is defined by 



U tj = (^2(i\H\d^(d^\H\j^J /{E d - Ei)(E d - Ej 



(1/5) 



Y^mid^id^H^/iEd-E,) 



J2{j\H\d v )(d v \H\j)/(E d -Ej 



(5) 



For the case described by Eq. (El), the interaction takes the form 
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Un = e( ri )e(rj) [P 2 (cos%) 2 - (1/5) j (6) 

where P2{@) = (3 cos 2 6 — l)/2 is the second-order Legendre polynomial. This, and the 
assumption that -Elfse is proportional to W, motivates the following choice for the functional 
form of -Elfse in terms of the local environment: 



-E-LFSE 

where 



n 1/2 



I>(r;)e(r,>(%; 



(7) 



M (fl) = [P 2 (cos%) 2 -(l/5)]. (8) 

Because the square root of the ligand-ligand sum is taken, this type of energy function is 
different from classical additive angular interaction potentials. I call it a "semiclassical" 
energy function, since the steps in its calculation are similar to those in the calculation of a 
classical energy function, but quantum mechanical effects are included in a systematic fash- 
ion. It applies to one spin component of a transition metal <i-shell; if both spin components 
contribute, then E'lfse is simply the sum of contributions from the two components. The 
form (0) is parallel in form to "many-atom" [13] and "embedded-atom" [FJ] energy func- 



tions, but these are not angle- dependent. Modifications of the embedded atom method fl5 
have included angular dependence, but without quantum-mechanical grounding, assuming 
angular forms very different from the present ones. 

In order to evaluate the accuracy of this functional form in the types of disordered local 
geometries that may be found in biomolecular environments, I have evaluated exact cluster 
energies (from the eigenvalues of Eq. (]!])) for an ensemble of transition- metal complexes 
having random bond lengths and angles. The transition metal ions have six neighboring 
ligands. The coupling strengths = e(r») in Eqs. (|]) and (0) vary randomly between 
and 2 (in arbitrary units), corresponding to distances varying from a short-range cutoff to 
infinity, and the orientations r, are chosen at random. In this way, a very broad range 
of environments, with effectively varying coordination numbers, is sampled. Semiclassical 
energy functions of the form fl7|), as well as classical energy functions, have been least- 
squares fitted to the exact d-electron energies of these clusters, for the ions Fe 2+ through 
Cu 2+ , taken in the high-spin configuration (Mn 2+ and Zn 2+ are not included, since their 
minority-spin ci-bands are empty and filled respectively, so -Elfse vanishes). In the fits, in 
addition to the ligand-ligand interaction terms, we include a constant term in the ligand- 
ligand interaction, as well as a sum of single-ligand terms. Figure 2a shows the fit for Cu 2+ 
obtained with the semiclassical energy function ([71). The energies are fit remarkably well, 
with the standard deviation of 0.16 being less than 10 percent of the typical values of I-ExfseI- 
Similar results are obtained for Ni. For Co 2+ , the fractional error is about 15 percent. For 
Fe 2+ , the magnitude of -Elfse is found to be an order of magnitude smaller than for Cu 2+ 
and Ni 2+ , and the fractional error resulting from using the potentials is about 50 percent; 
nevertheless, the absolute errors are about half of those for Cu 2+ and Ni 2+ . Figure 2b shows 



corresponding results for a classical potential for Cu + of the form y/eiy/ej sin 29, where the 



angular dependence is taken from recent simulations of cluster energetics [12] and the Je, 



dependence follows from dimensional analysis and the linear scaling of -Elfse with uniform 
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scaling the e^. We take this form to be typical of the treatment of transition metals in 
standard modeling packages in which simple plausible forms are assumed. The fit is much 
less accurate, with a standard deviation of 0.49. 

The energy function (^-|]) provides an optimal description of the G?-shell energetics in 
terms of two-ligand interaction interactions. To show this, I have fitted more elaborate 
potentials of the form (0) to the energy database, in which u(8) is represented by a sum 
of terms of the form cos n6, with n < 8. The results are shown in Fig. 3. The agreement 
between the optimized u{9) and the form (|8|) is almost exact for Cu 2+ and Ni 2+ , and very 
good for Co 2+ . For Fe 2+ , the absolute discrepancies are small, but the relative discrepancies 
are larger. Note that the shapes of u(8) as obtained here differ completely from the sin 2 28 
form of Ref. [ 12fl , which is shown by the dotted curve in frame (a). In addition, I have tried 
modified forms of Eq. (|7|), in which the square root is replaced by a power law dependence, 
so that an exponent of unity gives an additive potential. The minimum error is obtained 
with an exponent very close to 0.50, corresponding to the Eq. (|7]). Thus we have fairly 
definitively pinned down the functional form of the angular forces around these ions. We 
note that these results are also applicable to the low-spin versions of the ions, by simple 
addition of contributions from the two subbands. Then, for example, low-spin Ni 2+ becomes 
equivalent to high-spin Cu 2+ . 

The main chemical trend in u(6) with changing <i-count is a change in the magnitude 
of the potential, rather than its shape. The potentials for Ni 2+ and Cu 2+ are similar in 
magnitude, that for Co 2+ roughly a factor of two weaker, and that for Fe 2+ is weaker by 
an order of magnitude. The weakness of the Fe 2+ potential can be partly understood by 
analysis of the energetics of four-ligand complexes. For these, Hd, as given in Eq. ([!]), is 
a sum of four one- dimensional projection operators thus has rank four. One readily shows 
that all of its eigenvalues are nonnegative. This means that the lowest eigenvalue is zero, 
independent of the angular arrangement of the ligands. In the case of Fe 2+ , there is only 
one ci-electron, which resides in the orbital having the zero eigenvalue. Thus there are no 
angular interactions for Fe 2+ with four ligands. For cases with higher coordination, the 
lowest eigenvalue will still likely be close to zero unless the five contributing projection 
operators are orthogonal to each other. From the point of view of practical application, the 
variations seen in Fig. 3 suggest that the inclusion of angular forces for modeling Cu 2+ and 
Ni 2+ is crucial, but that the Fe 2+ ion (in high-spin configuration) might well be modeled 
with only radial interactions. 

These features can be used to explain the observed chemical trends in the relative sta- 
bility of square and tetrahedral structures in these systems. I have evaluated the energy 
difference AE between -Elfse between the square and tetrahedral coordination geometries. 
Comparisons between the exact values and those obtained by Eq. (|7|) and the empirical 
potential [|12|] are shown in Fig. 4, for the transition metal ions Fe 2+ through Cu 2+ . The 
empirical-potential results are much too small, but the basic trends of the exact results are 
also seen in the semiclassical results, with the square structure favored strongly for Ni 2+ 
and Cu 2+ . This trend is consistent with known structures of four-ligand transition metal 
complexes. Such complexes of Ni 2+ and Cu 2+ overwhelmingly adopt square coordination, 
in the absence of steric constraints, while Fe 2+ and Co 2+ generally have tetrahedral coor- 



dination |T^,|T8j. (We note that the experiments do not necessarily establish the sign of 



the electronic contribution AE calculated here for a given system, since direct electrostatic 
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interactions between the ligands tend to favor tetrahedral coordination; only the trend with 
varying ci-count is established.) The structural energies can be understood with the help of 
the potentials shown in Fig. 3. The minima at 0° and 180° favor the square structure in all 
cases, but are weaker for Fe 2+ and Co 2+ . In fact, the calculated values of AE correspond 
fairly closely to the strengths of the angular interactions. The energy differences are not, 
however, obtained quantitatively by the semiclassical energy function. The discrepancy lies 
mainly in the energy of the tetrahedral structure. For tetrahedral Co 2+ , for example, the 
semiclassical energy function underestimates |-Elfse| by about 20 percent. 

In summary, I have shown that a new semiclassical angular energy function, with a simple 
analytic angular dependence, describes the ligand-field stabilization energy for transition- 
metal ions in biomolecules remarkably well. The theoretical form for the angular dependence 
is strongly confirmed by analysis of a large computer-generated database of complexes. 
Analysis of the angular form of the interactions justifies the systematics of the relative 
stability of square and tetrahedral packing in terms of the behavior of the interactions at 0° 
and 180°. Incorporation of this form of energy function into existing bio molecular simulation 
packages should significantly enhance their reliability, and lead to new possibilities for design 
of metal-containing biomolecules. 
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FIGURES 



FIG. 1. Ligand-field splitting of Ni 2+ d-shell in square coordination. Only minority spin band, 
in high-spin configuration, is shown. 

FIG. 2. Accuracy test of semiclassical and empirical energy functions, in comparison with exact 
quantum-mechanical results for ligand-field Hamiltonian. Energy unit is average coupling of single 
ligand to transition-metal ci-shell. 

FIG. 3. Angular dependence of energy function. Solid lines: ten-parameter fit to exact energies. 
Dashed lines: derived angular function from Eq. (8). Function u{9) is dimensionless. Frame (a) 
Cu 2+ ; (b) Ni 2+ ; (c) Co 2+ ; (d) Fe 2+ . Dotted line in frame (a) is empirical energy function from 
Ref. [12], with magnitude adjusted for clear comparison. 

FIG. 4. Energy differences AE between square and tetrahedrally coordinated transition metal 
ions. Energy unit is coupling strength between single ligand and transition metal. Solid circles: 
exact treatment of ligand-field Hamiltonian. Open circles: semiclassical energy function (Eq. (7)). 
Triangles: empirical energy function (Ref. [12]). 
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